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SUMMARY 


OBJECTIVE 

Develop  a  spatial  smoothing  routine  that  accepts  arbitrary  ionospheric  models  as 
input.  Incorporate  the  smoothing  routine  into  a  sophisticated  ray  tracing  model  and 
state-of-the-art  ionospheric  model. 

RESULTS 

An  extension  of  a  smoothing  algorithm  previously  developed  by  A.  K.  Paul  of  the 
Naval  Ocean  Systems  Center  was  developed.  A  ray  tracing  model  developed  by 
R.  M.  Jones  and  J.  J.  Stephenson  was  implemented.  The  Utah  State  University 
Ionospheric  Model,  a  parameterized  version  obtained  from  D.  Anderson,  was  imple¬ 
mented. 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  iteration  process  required  to  obtain  the  coefficients  for  the  interpolation  uses 
too  much  computer  time;  a  faster  method  is  desirable.  The  ionospheric  model  needs 
to  be  changed  to  give  more  realistic  electron  densities  in  the  lower  altitudes;  this  will 
be  very  important  if  the  model  is  to  properly  account  for  absorption.  The  ray  trace 
model  has  some  difficulty  with  the  radical  changes  introduced  by  the  ionospheric 
model;  refinement  of  some  of  its  procedures  will  be  necessary. 


II 


INTRODUCTION 


The  Naval  Ocean  Systems  Center  (NOSC)  is  developing  the  High-Frequency  Bench¬ 
mark  Propagation  Analysis  Program,  or  HF  Benchmark.  This  program  employs  a 
sophisticated  ray  tracing  program  coupled  with  a  state-of-the-art  ionospheric  model. 

The  merging  of  these  two  major  components  depends  critically  on  the  development  of 
a  spatial  smoothing  routine  that  accepts  arbitrary  ionospheric  models  as  input.  Such  a 
spatial  smoothing  routine  has  been  developed  by  Paul  (1991).  This  report  describes  a 
computer  program  that  incorporates  an  extension  of  this  smoothing  algorithm  together 
with  the  Jones  and  Stephenson  (1975)  ray  tracing  routine  and  the  High  Latitude 
Ionospheric  Specification  Model  (HLISM).* 

Ionospheric  sounders  have  been  used  to  probe  the  ionosphere  for  more  than  40 
years.  Models  of  high-frequency  (HF)  propagation  are  based  on  the  data  acquired  by 
these  instruments.  But.  the  ionosphere  is  highly  variable.  Lacking  detailed  understand¬ 
ing  of  the  causes  of  ionospheric  variations  in  the  past,  and  being  limited  by  lack  of 
computing  power,  models  of  HF  propagation  have  relied  on  synoptic  data  to  develop 
global  maps  of  key  propagation  parameters  such  as  the  maximum  useable  frequency 
(MUF).  One  failing  of  this  approach  is  that  large  errors  in  prediction  of  critical  propa¬ 
gation  parameters  occur  under  many  important  conditions.  This  is  especially  true  at 
high  latitudes.  In  part,  errors  in  the  models  may  be  attributed  to  the  fact  that  they  are 
based  almost  exclusively  on  hourly  averages  taken  at  mostly  middle  latitude  stations. 
Another  important  contribution  to  errors  is  the  lack  of  reliable  models  of  intermittent 
phenomena  such  as  sporadic-E.  Recent  studies  of  data  measured  with  a  very  high  time 
resolution  show  that  short-term  fluctuations  of  the  ionosphere  may  be  large  and  that 
tilts  in  the  ionosphere  occur  often  (Paul,  1985). 

The  most  widely  used  HF  prediction  codes  in  use  today,  PROPHET  (Sailors,  1990) 
and  lONCAP  (Teters  et  al.,  1983),  are  largely  empirical,  based  on  semiempirical 
models  of  ionospheric  propagation.  These  programs  have  no  provision  for  first  princi¬ 
ples,  physics-based  calculations.  In  the  past,  first  principles  codes  were  too  slow,  com¬ 
plex,  and  user-unfriendly  to  address  the  general  class  of  problems  that  PROPHET  and 
lONCAP  attempted  to  address.  Today,  modern  desktop  computers  rival  the  mainframes 
of  the  past,  and  the  prospect  of  using  modern  first  principles  codes  in  routine  propaga¬ 
tion  problem  solving  is  a  reality.  Future  HF  prediction  models  will  be  hybrids  consist¬ 
ing  of  rigorous  physics-based  propagation,  expert  systems  to  monitor  and  interpret 
realtime  geophysical  and  ionospheric  sensors,  and  improved  man-machine  interfaces. 


*  S'Mircc  ?Vivatc  cc'nu7Uinicati<»n  vsiih  [).  And’'s«»n. 
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The  goal  of  the  work  described  in  this  report  is  to  develop  techniques  that  permit 
evaluation  of  different  propagation  models.  The  resulting  program  is  to  be  as  sophisti¬ 
cated  as  possible  with  concerns  for  computer  run  time  to  be  subordinate  to  the  use  of 
realistic  models  of  the  environment.  Employed  with  new  databases  of  high-latitude 
propagation  measurements,  this  model  will  be  used  to  evaluate  and  improve  the  faster- 
running  models  and  point  to  deficiencies  in  the  ionospheric  models. 

SPATIAL  SMOOTHING 

The  first  step  in  this  effort  was  to  develop  a  generalized  spatial  smoothing  routine 
that  will  accept  arbitrary  ionospheric  profiles.  A  candidate  spatial  smoothing  algorithm 
was  developed  by  Paul  (1991).  He  demonstrated  a  general  function  that  permitted  fit¬ 
ting  a  set  of  irregularly  spaced  discrete  data  with  satisfying  results.  Maps  of  foF2  over 
Europe  were  generated,  and  ray  tracing  in  the  subsequent  model  ionosphere  was  tested 
(using  no  magnetic  field  and  a  flat  earth).  This  algorithm  has  been  extended  to  use  a 
three-dimensional  array  of  ionospheric  profiles  of  electron  density  with  altitude  distrib¬ 
uted  in  a  grid  of  latitude  and  longitude. 

The  basic  procedure  in  this  algorithm  is  to  define  an  ionospheric  parameter  at  a 
series  of  coordinates.  Currently,  we  chose  this  parameter  to  be  the  logarithm  of  the 
electron  density.  A  single  data  point  is  characterized  by  a  function,  /„  where 


/,(.v,y,z)  =  1  + 


(^1  •  ('?)■  •  fe)' ' 


(1) 


1  + 


and  the  vrlue  at  a  point  {.x,y,z)  due  to  all  of  the  data  points  in  the  set  is  given  by 


p{x,y,z)  =p„  +  p,  nf,{x,y\z). 


(2) 


The  quantities  A;-  A/,  and  r^,  are  called  the  control  distances  and  basically  determine 
the  range  of  influence  of  individual  data  points.  We  have  found  that  these  control  dis¬ 
tances  should  be  about  80  percent  of  the  separation  of  the  data  points  in  each  direc¬ 
tion.  The  quantities  p„  and  p^  are  scaling  parameters.  We  use  them  to  normalize/!. 
This  improves  the  convergence  of  the  iterations  used  to  find  the  values  of  the  unknown 
coefficients,  a,. 

Modifications  fo  the  original  procedure  were  required  to  consistently  achieve  conver¬ 
gence  of  the  coefficients.  Given  that  there  is  no  functional  form  to  impose  on  the  itera¬ 
tions.  the  following  procedure  is  used.  The  minimum  value  of  the  electron  density  data 


is  determined  (po)  and  subtracted  from  all  of  the  values.  These  data  are  then  normal¬ 
ized  using  the  peak  value  of  the  offset  data  (p,).  Equation  2  is  used  to  compute  a 
value  at  each  of  the  input  locations.  A  correction  term  is  defined  equal  to  the  differ¬ 
ence  between  the  computed  and  input  values  of  the  logarithm  of  the  electron  density  at 
the  input  locations.  If  the  largest  of  these  correction  terms  is  less  than  a  specified 
amount,  then  the  iterations  are  terminated.  Otherwise,  the  coefficients  are  adjusted  by 
some  fraction,  and  the  calculations  are  repeated.  Currently,  we  are  using  a  maximum 
difference  between  these  calculations  and  the  input  data  of  0.1,  which  amounts  to  10 
percent  of  the  logarithm  of  the  peak  electron  density  and  seems  to  give  reasonable 
results  in  about  30  iterations.  We  have  run  the  algorithm  using  0.01  or  1  percent  of 
the  peak  value  of  the  logarithm.  This  took  120  iterations  and  produced  no  discernable 
difference  in  the  ray  paths. 

The  large  number  of  data  points  required  to  give  a  reasonable  representation  of  the 
ionospheric  variation  along  a  propagation  path  seems  to  require  that  the  fraction  used 
to  adjust  the  coefficients  be  quite  small  (0.05  in  the  current  configuration  of  the  pro¬ 
gram).  This  small  value  prevents  oscillatory  behavior  of  the  iterative  corrections,  but 
requires  a  lot  of  iterations.  We  have  implemented  two  procedures  that  improve  the 
efficiency  of  the  iterations.  After  much  experimentation,  we  found  that  the  initial 
estimate  for  the  u,  should  be  some  fraction  of  the  quantity  (/,-!.  Currently  we  use  0.6 
for  this  fraction.  The  time  required  to  calculate  equation  2  is  reduced  by  defining  a 
factor  that  controls  the  range  of  influence  of  each  data  point.  It  can  be  seen  in  equa¬ 
tion  2  that  this  is  useful  because  the  influence  of  each  input  data  value  decreases  as 
the  square  of  the  distance.  Data  points  outside  of  the  range  of  influence  of  an  interpo¬ 
lated  point  are  not  included  in  the  calculation.  The  improvement  in  computation  time 
is  achieved  by  eliminating  the  large  number  of  floating-point  divisions  found  in  equa¬ 
tion  2  and  avoiding  calculation  of  terms  that  will  not  contribute  to  the  total.  In  the  cur¬ 
rent  program,  this  range  of  influence  is  six  times  each  of  the  control  distances.  This 
procedure  also  reduces  the  time  required  to  generate  interpolated  values. 

Three  subroutines  are  used  to  implement  the  interpolation.  The  routine  named 
AKP_MODEL  takes  the  initial  path  parameters  and  sets  up  the  array  of  data  points  by 
calculating  the  locations  of  the  ionospheric  profiles.  The  routine  AKP  COEFF  finds  the 
quantities,  p^,  and  p,.  and  performs  the  iterations  to  determine  the  coefficients,  o,.  The 
routine  AKP  XYTHPH  does  coordinate  conversion  between  the  magnetic  dipole  coordi¬ 
nate  system  used  in  the  ray  trace  routines  and  the  quasi-Cartesian  coordinate  system 
used  in  the  interpolation  routine.  This  routine  also  calculates  the  derivatives  of  the 
ionospheric  parameters  due  to  the  spatial  variation  of  the  ionosphere.  The  equations 
for  performing  these  calculations  are  summarized  in  appendix  A. 
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IONOSPHERIC  MODEL 


The  ionospheric  model  currently  used  in  the  program  is  the  Utah  State  University 
(USU)  ionospheric  model,  a  parameterized  version  which  we  obtained  from  Anderson.* 
This  model  is  called  the  High  Latitude  Ionospheric  Specification  Model  (HLISM).  We 
have  developed  a  subroutine  driver  for  the  model  called  USU_MODEL.  This  routine  is 
called  by  AKP_MODEL  to  generate  ionospheric  profiles  at  points  500  km  apart  along 
the  propagation  path  (defined  by  the  transmitter  and  the  receiver''  and  at  points  500 
km  away  on  either  side  of  the  path  along  a  direction  perpendicular  to  the  propagation 
path. 

We  chose  the  HLISM  because  it  is  considered  one  of  the  best  high-latitude  models 
a\ailable.  Unfortunately,  it  is  a  model  for  positive  ion  densities  in  the  altitude  range 
from  100  to  800  km.  We  invoke  charge  neutrality  and  set  the  electron  density  equal  to 
the  total  ion  density.  This  is  valid  only  where  the  negative  ions  can  be  ignored.  The 
model  tends  to  generate  (correctly)  large  ion  densities  at  the  bottom  of  the  ionosphere 
that  would  normally  be  neutralized  by  the  presence  of  negative  ions  as  well  as  elec¬ 
trons  at  those  altitudes.  For  the  moment,  we  have  made  an  adjustment  to  the  model  as 
follows.  We  calculate  the  electron  density  down  to  160  km  and  then  append  a  value  of 
10^  at  100  km  and  lO""  at  60  km  to  the  bottom  of  the  profile.  Under  nighttime  condi¬ 
tions,  these  values  may  cause  a  slight  kink  in  the  electron  density  profile.  The  adjust¬ 
ment  is  required  to  provide  a  gradual  entry  into  the  ionosphere  for  the  ray  tracing  rou¬ 
tine.  We  found  that  the  unrealistically  high  values  of  electron  densities  at  the  bottom 
of  the  ionosphere  were  causing  numerical  problems  in  the  ray  tracing.  We  plan  to 
obtain  and  implement  an  electron  density  model  either  as  a  modification  of  the  USU 
model  or  as  a  replacement  for  it. 


RAY  TRACING 

We  have  implemented  the  Jones  and  Stephenson  (1975)  ray  tracing  model.  This  is 
a  versatile  program  with  full  allowance  for  externally  specified  models  of  the  electron 
density,  collision  frequency,  and  geomagnetic  field.  We  have  made  every  effort  to 
retain  this  flexibility  in  our  implementation.  In  particular,  the  input  requirements  are 
nearly  identical  to  the  original  Jones  and  Stephenson  version  with  our  model  being 
called  out  as  “AKPMODEL”  in  the  specification  of  the  ionospheric  model.  The  output 
from  the  original  model  is  retained.  We  have  added  additional  parameters  for  input. 
These  are  summarized  in  table  1 . 


*  Private  communication  with  D.  Anderson,  19'>1. 


Minor  procedural  modifications  have  been  made  to  facilitate  running  this  program 
on  a  PC,  but  the  FORTRAN  used  in  the  program  is  not  specific  to  that  platform.  In 
fact,  an  identical  capability  exists  for  a  VAX. 


Table  1.  Ray  traee  inputs  for  AKP  model. 


— 

9 

This  value  used  to  signal  that  the  a,  have  already  been 
computed. 

111 

1 

Number  of  the  month 

112 

31 

Day  of  the  month 

113 

1991 

Year 

114 

1200 

UT  in  HHMM  format 

115 

100 

10  cm  flux 

116 

1 

117 

1 

Sign  of  By;  1  for  -1  for 

118 

160 

Minimum  USU_MODEL  height  to  use;  km 

119 

600 

Maximum  USU_MODEL  height  to  use;  km 

120 

0.8 

Horizontal  control  distance  factor 

121 

0.8 

Vertical  control  distance  factor 

■B 

0.6 

Initial  scale  factor  for  Oi 

123 

6.0 

Range  of  influence  factor 

124 

0.0625 

Iteration  adjustment  factor 

125 

0.1 

Maximum  error  to  terminate  iterations 

126 


51 


Maximum  number  of  iterations 


SAMPLE  PROBLEM 


A  sample  problem  is  presented  here  to  illustrate  the  unification  of  the  individual 
models.  The  path  is  defined  in  geographic  coordinates  from  40°N,  90°W  to  80°N, 
90‘’W.  The  geographic  bearing  angle  is  0°  an<l  the  path  length  is  4400  km.  The 
ionospheric  profiles  are  generated  for  15  July  1991  at  1800  hours  UT.  The  10-cm  flu.x 
is  200,  the  kp  is  5,  and  the  direction  of  the  sun’s  magnetic  field,  By,  is  positive.  Con¬ 
tour  plots  of  the  logarithm  of  the  electron  density  in  vertical  planes  along  geographic 
bearings  of  -5°,  0°,  and  5°  are  shown  in  figure  1.  Clearly,  the  ionosphere  is  changing 
in  all  three  dimensions.  In  particular,  we  note  a  broad  ma.ximum  near  2  Mm,  followed 
by  a  constriction  of  the  contour  lines  near  3.2  Mm.  Figure  2  shows  the  resulting  ray 
paths.  The  top  panel  (figure  2a)  shows  the  projection  of  the  ray  paths  onto  the  verti¬ 
cal  plane  passing  through  the  transmitter  and  the  receiver.  The  rays  with  small  circles 
on  their  ends  are  rays  that  penetrate  the  ionosphere,  and  we  see  that  this  penetration 
occurs  where  we  see  a  reduction  in  electron  density  in  figure  1.  The  bottom  panel  of 
figure  2  (figure  2b)  shows  the  projection  of  the  ray  paths  onto  the  ground  with  loca¬ 
tions  above  the  horizontal  axis  being  to  the  west  of  the  propagation  path.  We  see  that 
the  east-west  variation  of  the  ionosphere  causes  some  deflection  of  the  ray  paths. 


CONCLUSION  AND  RECOMMENDATIONS 

The  iteration  process  required  to  obtain  the  coefficients  for  the  interpolation  model 
requires  excessive  computer  time,  and  the  interpolation  is  quite  slow.  A  faster  method 
is  highly  desirable.  The  current  interpolation  model  still  requires  too  much  manual 
adjustment  to  be  routinely  used.  The  ionospheric  model  needs  to  be  changed  to  gi\c 
more  realistic  electron  densities  in  the  lower  altitudes.  This  wilt  be  very  important  if  the 
model  is  to  properly  account  for  absorption.  The  ray  trace  model  has  some  difficulty  with 
the  radical  changes  introduced  by  the  ionospheric  model,  and  refinement  of  some  of 
its  procedures  will  be  necessary. 
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APPENDIX  A:  COORDINATE  TRANSFORMATIONS. 


Coordinate  conversion  is  needed  because  the  ray  trace  program  works  in  spherical 
coordinates  and  the  interpolation  routine  uses  Cartesian  coordinates.  A  position,  P,  on 
a  sphere  defined  by  a  colat'tude,  0  and  longitude,  (p  must  be  converted  to  a  distance, 
A',  along  a  reference  great  circle  path  from  the  transmitter  and  a  distance,  y,  perpen¬ 
dicular  to  that  great  circle.  We  chose  this  reference  great  circle  to  be  the  one  which 
passes  through  both  the  transmitter  and  the  receiver.  The  relevant  geometry  is  shown 
in  figure  A-1.  The  transmitter  position  is  labeled  T  and  is  located  at  colatitude  0,  and 
longitude  0,.  The  great  circle  between  the  transmitter  and  the  receiver  (not  shown) 
makes  an  angle,  with  respect  to  north.  The  position  for  w'hich  we  need  interpolation 
coordinates,  a  and  y,  is  labeled  P.  In  the  figure,  50  is  merely  the  difference  in  longi¬ 
tude  between  the  transmitter  and  P.  Note  that  the  line  segment,  y,  is  perpendicular  to 
the  reference  great  circle. 


Figure  A-1.  Illustration  of  spherical  coordinates  to  interpolation 
coordinates. 


A-1 


The  law  of  cosines  for  sides  gives  the  great  circle  distance,  r,  and  the  bearing 
angle,  a ,  from  th'^  transmitter  to  P\ 


cosr=  cos  6  COS0,  +  sin0  sin  9,  cos  6cp 


cos  6  -COS0,  COST 

cosa  =  - r— ^ ^ - 

sin  6^,  sinr 


The  law  of  sines  gives 

s\nO  sin(50 

sin  O!  = - — 

sin  r 


Referring  to  the  right  triangle  defined  by  the  sides  labeled  r,  x,  and  y,  the  law'  of 
cosines  for  sides  gives  us 


cos  r 

cos  X  -  — — . 
cos  V 


The  law  of  sines  gives 


siny  =  cos(/i  -  o; )  cosy 

and  the  law  of  cosines  for  angles  gives  us 

cos  (6~Oi) 

siny  =  — — c - 

cosy 

Finally,  the  law  of  sines  gives 
sinA'  =  siny  sin r. 


These  equations  can  be  combined  to  eliminate  ot,  r,  and  y ; 

siny  =  cos 0  sin0,  sin/3  -  sin0  (cos/3  sin 50  +  cos0t  sin/3  cos 50) 


sinA  = 


cosf^  sin6>i  cos/3  +  sin0  (sin/3  sin 50  -  cos0i  cos/3  cos 50) 

cosy 


cos  A 


cosO  cosOi  +  sin0  sin0,  cos  50 
cosy 


A-2 


The  ray  tracing  routine  also  requires  the  derivatives  of  the  ionospheric  parameters 
with  respect  to  the  spatial  dimensions.  The  necessary  equations  are  listed  below. 


=  -s\n0  sinO,  sin/S  -  cos  6^  (cos sin  (30  +  cos^,  sin/3  cos  50) 
=  -  sin  0  (cos/3  cos 50  -  cos0,  sin/3  sin  50) 


-sin0  sin^,  cos/3  +  cos 6  (sin  ^  sin 50  -  cos^t  cos/3  cos 50) 


dO 

d  sin.r 

sin/3(sin/3  sin  50 

cosy 

+  cos^,  cos/3  sin 50) 

(30 

cosy 

b  cos  X 

-sin0  COS0,  +  COS0  sin  6, 

cos  50 

b0 

cosy 

b  COS  .r 

sin0  sin0,  sin  50 

b(p 

cosy 

d  sinx 

b  sin  x 
_  1 

sinA  siny  c/siny 

d0 

b0 

cos^y 

d0 

d  sin.v 

b  sinx 

sin  A  siny  c/siny 

c/0 

1 

50 

cos^y 

d(P 

d  cosx 

b  COS  a: 

—  X 

cos  A  siny  c/siny 

d0 

b0 

cos^y 

d0 

d  cos  X 

5  cos  A 

cos  A  siny  d  siny 

d(p 

_  ■+■  ' 
50 

cos^y 

d(p 

d  siny 
dO 

d  siny 
£/0 


5  sin.t 
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